План митапа

Типы пространственных данных;
Особенности геоданных;
R и пространственный анализ;
Практика - находим геоданные разных форматов и учимся с ними работать;
Что можно почитать и посмотреть по теме.

1. Типы пространственных данных

а. Векторные данные

Геометрии состоят из точек. Конвенционально x = широта, y = долгота. Кроме этого могут использоваться z = высота, M = мера изменчивости точки;
Виды геометрий:

+++ GEOMETRYCOLLECTION +++

Форматы записи данных:

Well-Known Text (WKT) - Запись точек вектором из геометрий

Well-Known Binary (WKB) - Запись координат в бинарных значениях. Используется в базах данных (т.к. увеличивает скорость работы с данными), но нечитаем и непонятен для человека.

Посмотрим в первом приближении:

library(sf)

library(rnaturalearth)
world <- ne_countries(returnclass = "sf")
world <- world[1:5,c(4, ncol(world))]

# Для иллюстрации
world_1 <- world
world_1$geometry <- as.character(world_1$geometry)
as.data.table(world_1)

b. Растровые данные

Матрица значений пикселей пространственной области. В растрах хронят космические снимки и базы различных геологических данных. Также как и вектор - предназначены для определённого круга задач;
Основной формат: .tif
Основные библиотеки для обработки в R: raster и stars

Сравним:

На этом и последующих занятиях мы будем говорить только о векторных данных.

Главное, нужно помнить, что векторные пространственные данные - это (чаще всего) просто точки на плоскости! Создадим собственные пространственные данные:

p1 = st_point(c(7,52)) # Рисуем точку 1
p2 = st_point(c(-30,20)) # Рисуем точку 2
sp_obj = st_sfc(p1, p2, crs = 4326) # Преобразуем в пространственные данные
plot(sp_obj)

# Вы восхитительны!

…

2. Особенности геоданных

Вернёмся к векторным геометриям:

world

Coordinate Reference Systems (CRS) / Пространственная привязка и её проекции

Что это такое?

Для иллюстрации обратимся к прекрасной работе Тасс.

Картинки чтобы окончательно закрепить идею проекции.

Как посмотреть проекцию в CRS?

st_crs(world)

## Coordinate Reference System:
##   User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

На самом деле всё не так сложно.

Пространственный объект можно перевести в другую CRS сразу и быстро, задав через st_transform() нужную проекцию одним из четырёх форматов (просто гуглите то, что нужно для вашего пространства):

Код EPSG (например, “4326”)
Формула PROJ4 (например, “+proj=longlat +datum=WGS84 +no_defs”)
Строчка WKT
Просто переносом проекции с другого пространственного объекта ( df1 <- st_transform(df1, st_crs(df2)) )

Например:

Afg = world[1,]
plot(Afg)

# Change crs
Afg <- st_transform(Afg, 2264)
plot(Afg)

# Change crs
Afg <- st_transform(Afg, "+proj=longlat +datum=WGS84 +no_defs")
plot(Afg)

В геоданных есть ещё множество подводных камней, но этот - ключевой и его необходимо запомнить.

3. R и пространственный анализ

Почему R?

Библиотеки для работы с геоданными в R

(несть им числа, но основные)

sf - state-of-the-art пространственного анализа в R.

Немного истории:

4. Практика

Что мы будем делать - делать карту мороженного в Израиле! :)

Я исхожу из того, что вы знаете dplyr или data.table на базовом уровне, но если у вас вызывает затруднения именно манипуляции с данными - обязательно задавайте вопросы!

Мы отработаем работу с данными на трёх “классических” форматах для векторных данных: csv, geojson, shape-file.

A. Датасет с мороженками в csv

Загрузим нужные библиотеки:

library(sf)
library(geojsonsf)
library(geojsonio)
library(readr)
library(data.table)

Загрузим данные. К сожалению, в R есть известная проблема с кодировками и чтобы долго не решать проблему с крокозябрами - уберём колонки названий точек мороженного на иврите:

isr_icecream <- fread('https://raw.githubusercontent.com/AmitLevinson/Datasets/master/golda/golda_locations.csv', encoding = 'UTF-8')
isr_icecream[,1:2 := NULL]
isr_icecream

Превратим в sf тип данных:

isr_icecream <- st_as_sf(isr_icecream, dim = "XY", remove = T, na.fail = F, coords = c("lon", "lat"), crs = "+proj=longlat +datum=WGS84 +no_defs")
isr_icecream

plot(isr_icecream)

B. Границы Израиля в Shape-file

Идём сюда, ищем Израиль и загружаем архив. Что мы видим внутри?

.shp – главноеый файл с геометриями

.dbf

.shx

.prj

(но все они важны!)

По этой причине, когда вам нужно переслать shape-file отправляйте архив из всех четырёх файлов.

Загрузим его и отобразим:

isr_border <- st_read('/cloud/project/ISR_adm/ISR_adm0.shp')

## Reading layer `ISR_adm0' from data source `/cloud/project/ISR_adm/ISR_adm0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 34.26801 ymin: 29.49708 xmax: 35.90094 ymax: 33.36403
## CRS:            4326

isr_border

plot(st_geometry(isr_border))

C. Точки городов Израиля в GeoJson

Json с геоданными. Достаточно тяжёлый с точки зрения размещения данных, но часто встречающийся формат.

И здесь стоит отвлечься и ознакомиться с “Википедией”" от мира пространственных данных - Open Street Map.

Если это “Википедия”, то с неё можно выкачать данные? Да. По тэгам…

…

Тэги городов.

Скачаем GeoJson с городами. Загрузим его. Внутри будет огромное количество колонок, которые сейчас нам не нужны. Оставим первую и последнюю:

isr_city <- st_read('/cloud/project/city.geojson')

## Reading layer `city' from data source `/cloud/project/city.geojson' using driver `GeoJSON'
## Simple feature collection with 51 features and 225 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 33.80459 ymin: 30.61197 xmax: 36.56875 ymax: 32.92817
## CRS:            4326

isr_city <-  isr_city[,c(1,ncol(isr_city))] 
isr_city

plot(isr_city)

Альтернативный способ через пакет osmdata (спасибо Philipp Pyshny за наводку)

library(osmdata)

admin_osm = opq(bbox = 'Saint-Petersburg') %>%
  add_osm_feature(key = "admin_level", value = '5') %>%
  osmdata_sf()

# Но есть проблемы с кодировкой
iconv(admin_osm$nodes$tags, from="UTF-8", to="UTF-8") 
iconv(admin_osm$nodes$tags, from="UTF-8", to="cp1251")

admin_osm[["osm_multipolygons"]]

Объединяем слои и строим свою первую карту!

# Устанавливаем одинаковые проекции для всех объектов: 
isr_border <- st_transform(isr_border, 4326)
isr_city = st_transform(isr_city, st_crs(isr_border))
isr_icecream = st_transform(isr_icecream, st_crs(isr_border))

# Делаем карту (st_geometry() - вывести чистую геометрию, без других колонок или значений в объекте)
{plot(isr_border %>% st_geometry())
plot(isr_city, col = 'red', add = T)
plot(isr_icecream, col = 'blue', add = T)}

Ура :)

Теперь сохраним данные (здесь есть свои подводные камни):

# write_sf(isr_border, "isr_border.shp", append = T, layer_options = "ENCODING=UTF-8")

Загрузим в качестве проверки и убедимся, что всё работает

# double_isr_border <- st_read("isr_border.shp")
# plot(st_geometry(double_isr_border))

Бонус

Построим интерактивный график с помощью очень простой по своему синтаксису библиотеки mapview. Туториал по ней можно найти здесь. mapview предназначен для карт небольшого размера.

ЗАПУСКАТЬ ОСТОРОЖНО!!!

(на Windows возможны вылеты RStudio)

library(mapview)
mapview(isr_border) + 
  mapview(isr_city %>% st_geometry(), col.regions = 'red', legend = FALSE) +
  mapview(isr_icecream %>% st_geometry(), col.regions = 'blue', legend = FALSE)